
import numpy as np
import matplotlib.pyplot as plt
# [ y_1, y_2 ]^\top
def f(x): # x 是一个 2 x 1 的列矩阵
    return np.mat([[x[0,0] * np.cos(x[1,0])-4], # 矩阵定义： np.mat
            [x[1,0] * x[0,0] - x[1,0] - 5]]) # 矩阵元素：x[i,j]

# partial derivatives: Jacobi matrix
def df(x):
    return np.mat([[np.cos(x[1,0]), -x[0,0] * np.sin(x[1,0])],[x[1,0], x[0,0]-1]])

def newton2d(x, f, precision):
    mDeltax = 10 * precision # module of Delta x
    iter = 0
    while( mDeltax > precision ):
        Deltax = - df(x).I * f(x) # Delta x = - J^{-1} * y
        mDeltax = np.sqrt((Deltax.T * Deltax)[0,0]) # module of Delta x
        #print("mDeltax = ", mDeltax)
        x = x + Deltax
        iter += 1
        #print("new x = ", x.T, " new f = ", f(x).T)
    return x, f(x), iter

x = np.mat([[1],[1]])
z = newton2d( x, f, 0.01)
print("After ", z[2], "iterations, root = ", z[0].T, ", f(root) = ", z[1].T)